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Abstract 

This manual describes the basic objectives, functionalities and uses of 
the toolbox for Maple (Maplesoft™ ) called Quantavo. It is intended to 
facilitate calculations both symbolically and numerically related to Quan- 
tum Optics. In particular the evolution, measurement and entanglement 
properties of quantum states in the Fock basis can be simulated with it. 
It is provided to the community as a free open source module. 
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1 Objectives 

The toolbox is made to be used with the following approach: 

Declare an initial state, let it evolve through a quantum optical circuit involving 
linear optics (LO) and measurements and finally ask various questions about 
the structure, entanglement, and properties of the final state. 



The toolbox intends to do so providing: 

A) A framework in which to declare, manipulate and characterize quantum 
states of light (finite number of modes, and finite dimensional). 

B) Procedures that implement linear operations or whole linear optics circuits 
on our states such as: 

• Beam Splitters (BS) 

• Phase Shifters (PS). 

• arbitrary unitary transformations of the modes. 
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C) Procedures that implement arbitrary measurements (both projective or gen- 
eralized positive operator valued measures (POVM)). 

D) Procedures to determine probabilities and expected values for projective 
measurements and POVM measurements. 

E) Procedures to trace out measured or inaccessible modes. 

F) Procedures to calculate different entanglement measures such as: 

• Entropy of entanglement 

• partial trace, norm — > Negativity, Logarithmic Negativity 

G) Access to properties such as the mean Energy of given states. 

H) Extract and display lists of coefficients along with their indexes like: RTy/2 
|001)(001| 

H) Tools to easily plot states and density matrices. 

2 Getting Started 

Some previous knowledge about Maple from Maplesoft™ is required to use this 
tool. However, the brief "Take a Tour of Maple" should suffice to get started. 
This toolbox works with Maple 9.5, Maple 10 and Maple 11. 

2.1 Definitions and Notation 

Procedures: Formally, in Maple, a procedure definition is a valid expression 
that can be assigned to a name. The procedures we will use can be thought of 
as a set of "rules" that generally receive one or more inputs and return one or 
more outputs. We will write them in typewriter face. For example the proce- 
dure IsHermitian(M) evaluates if a Matrix M describing a density matrix is 
hermitian. A dictionary with all the procedures involved in Quantavo can be 
found in the appendix. 

Modules: Modules are repositories of procedures. By loading a module we 
can use its procedures in the Maple worksheet. For example Quantavo is the 
module containing the procedures we will use. 

objects: The objects in which we will encode our quantum states will be writ- 
ten in bold face. These will include vec, mat, matcol and poly and will be 
introduced later. 

d and K: Throughout the manual, 'if' will stand for the number of modes 
and 'cf for the dimension of each mode. If considering the photon number de- 
gree of freedom, then 'd — V will be the maximum number of photons in any 
mode. It is important to keep track of the value of these two global variables 
throughout the worksheet as they play an important role in the labelling of 
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the optical modes and translation procedures. These variables can always be 
updated and displayed with the procedure: f indKnd(State). 

2.2 Loading the Modules 

All necessary files can be found at, 

http : //www. imperial . ac .uk/quantuminf ormat ion/research/downloads 

1. uncompress the file Quantavo.zip. 

2. Save the folder QUANTAVO to a given directory. It should contain the 
files Quantavo.mpl and Quantavo_Example_Worksheet.mw. 

3. One can start opening the worksheet Quantavo_Example_Worksheet.mw. 

4. To use the module in a new worksheet execute the following commands 0: 



> with (Linear Algebra): 

> read "Quantavo.mpl" ; 

> with(Quantavo); 



It should return a list of all the procedures available 



[APD, BS, BuildUnitary, CoherentState, DP, Dbra, Dbraket, DeltaK, 
Dket, Dstate, Energy, Entropy, EvalState, IdentityState, IsHermitian, 
IsNormalized, LogNegativity, Negativity, POVMresult, PS, PlotState, 
Probability, Project, SqueezedVac, StateApprox, StateComplexConju- 
gate, StateMultiply, StateNorm, StateNormalize, StatePartialTranspose, 
StateSort, StateTrace, TensorProduct, TensorVac, Traceout, Trim, 
UnitaryEvolution, Vac, findKnd, indexstate, mat2matcol, mat2poly, 
matcol2mat, matcol2poly, modesmatcol, myBS, poly2matcol, poly2vec, 
vec2mat, vec2matcol, vec2poly] 



5. You are ready to use Quantavo ! 



Note that You may also save the module to your Maple library. To do so, 
visit the maple help on module, savelibname and savelib. 



1 note that if the new file is in a different directory it should be, 
>read "/path-to-folder/QUANTAVO/Quantavo.mpl" ; 
or under MS Windows, 

>read "C:\\path-to-folder\\QUANTAVO\\Quantavo.mpl" ; 
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3 Toolbox 



To run Quantavo one has to load the following modules: 

• Linear Algebra (Linear Algebra package from Maplcsoft™ built in Maple) 

• Quantavo (General toolbox) 

Additionally, if we want to plot our states we will also need the module 

• geom3d (geometry package from Maplesoft™ built in Maple) 

• plots (plotting package from Maplesoft™ built in Maple) 

3.1 Objects and Operations 

When quantum optical states have a few modes and live in high dimensions, the 
matrices or vectors describing them very soon become intractable. To mitigate 
this difficulty, the procedures from Quantavo store and manipulate only the 
non-zero elements in the description of our states. Our main objects will be 2 
column and 3 column matrices. 

3.1.1 Pure States 

Two column objects will describe vectors in Hilbert space of the form: 

d-l 

\<t>) = f(ni,n 2 ,...,n K )\n 1 ,n 2 ,...,n K ) 

ni,n2,...,n K 

Where all indices "m" range from — > (d — 1). Pure state vectors will be 
encoded in 2 column matrices that we will call trimmed vectors containing only 
non-zero entries. They will only contain non-zero entries. These objects will be 
named as a short vector: "vec" and will have the following appearance: 



1 


[0,0,0] 


A 


[1,1,0] 


A 2 


[2,2,0] 


A 3 


[3,3,0] 



For each row, the second column will be a list with the number of photons in 
each mode, This way [0,1,2] stands for 1 012) . the first column will contain the 
coefficient associated with this ket. The whole will describe the linear supper- 
position of all these kets with their coefficients, therefore the above -0 describes 
the unnormalized quantum state: 

3 

M = £A>,n>. 

n=0 
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3.1.2 Mixed States 



Two objects will be used to display density matrices. The first one is a square 
matrix, with as little zero entries as possible. This matrix will be called "mat" 
and will have for the above state the following form: 






[0,0,0] 


[1,1,0] 


[2,2,0] 


[3,3,0] 


[0,0,0] 


1 


A 


A 2 


A 3 


[1,1,0] 


A 


AA 


A A 2 


A A 3 


[2,2,0] 


A 2 


A 2 A 


A 2 A 2 


A 2 A 3 


[3,3,0] 


A 3 


A 3 A 


A 3 A 2 


A 3 A 3 



Another object that can also describe a density matrix will be a 3 column 
matrix, or trimmed density matrix. It will be named as a short column matrix: 
"matcol". The second column will be a list with the number of photons in 
each mode of the ket. This way, [0,1] in the 2nd column describes |01). The 
third column will be a list with the number of photons in each mode of the bra; 
therefore [2,1] in the 3rd column describes (21|. Finally The first column will 
have the coefficient associated with this |ket)(bra| . The whole will describe 
the non zero elements of the density matrix. As an example consider the above 
state which will be: 



Pmatcol 



1 


[0 


0,0] 


[0,0,0] 


A 


[0 


0,0] 


[1,1,0] 


A 2 


[0 


0,0] 


[2,2,0] 


A 3 


[0 


0,0] 


[3,3,0] 


A 


[1 


1,0] 


[0,0,0] 


A A 


[1 


1,0] 


[1,1,0] 


A A 2 


[1 


1,0] 


[2,2,0] 


A A 3 


[1 


1,0] 


[3,3,0] 


A 2 


[2 


2,0] 


[0,0,0] 


A 2 A 


[2 


2,0] 


[1,1,0] 


A 2 A 2 


[2 


2,0] 


[2,2,0] 


A 2 A 3 


[2 


2,0] 


[3,3,0] 


A 3 


[3 


3,0] 


[0,0,0] 


A 3 A 


[3 


3,0] 


[1,1,0] 


A 3 A 2 


[3 


3,0] 


[2,2,0] 


A 3 A 3 


[3 


3,0] 


[3,3,0] 



Summarising, we will use mainly 3 objects: vec, mat and matcol. 
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A fourth object less commonly used is a polynomial representation of the state. 
It is a polynomial in the mode operators that define the state. A general state 
would then be: 

d-i 

P = f{n 1 ,n 2 ,...,n K ,m 1 ,...,m K ) a\ ni a\ n \..a^ K |0)(0| b^b™\..V 

ni,n 2 , ...,n K 
mi,TO 2 , ...,m K 

and its description as a poly object: 

poly(p)= Yl f(n 1 ,n 2 ,...,n K ,m 1 ,...,m K ) a n ^a n 2 \..a n K K b^b™\..b™> 

ni,n 2 , ...,n K 
m 1 ,m 2 , ...,m K 

(note that the commutation relations are not taken care of in poly objects). 

3.2 Declare, Propagate, Measure and Ask 

Quantavo contains various procedures that allow problems to be formulated 
as follows: 

1. Declare the initial state. 

2. Apply different transformations to it (Beam Splitter, Phase Shifter, Arbi- 
trary Unitary, ...). 

3. Measure certain modes (and trace out the inaccessible ones), find out 
probabilities. 

4. Ask different questions about the properties of the state: Display, Plot, 
evaluate certain measures of Entanglement, etc.. 

Additionally, the order in which we use these procedures can be changed 
and the questions in item (4) can be formulated at any intermediate time. In 
addition, more states can be tensored or added at later times. Finally there 
are procedures to interconvert mat, vec, matcol and poly. Let us then give a 
more detailed description of these four basic steps. 

3.2.1 Declaration 

There are different ways to declare a state, all depending on its characteristics. 
Pure States 

For pure states of known functional form. That is, if we have a state of the 
form: 

d 

|V>) = Y f{ni,n 2 ,...,n K )\n 1 ,n 2 ,...,n K ) (1) 

n 1 ,n 2 ,...,n K 
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and we know explicitly /(ni,ri2, we may use the following structure: 

Declare the number of modes "if" , the maximum number of photons "c?" , and 
make a loop to declare the elements. For example: 



d:= 4; 
K := 3; 

V:=Matrix(cF,2): 



There are now 2d elements to be specified. 



for i from 1 to d do 
for j from 1 to d do 
for k from 1 to d do 

V[i,l]:=/(i,j,k): 
V[i,2] :=[i,j,k]: 

end do: 
end do: 
end do: 



Where / is the function in eq.JT]). We may use the procedure deltaK(ij) if 
we need a Kronecker delta in our definition. Executing the above loop will de- 
clare a matrix V of size d K x 2 that has hopefully many zero entries. To get rid 
of the zero entries and convert this object into a vec object (or trimmed vector) 
we will use the procedure Trim: 



Vl:=Trim(V): 



Another way to proceed is to declare the object vec directly with an appro- 
priate function. Consider as an example the state: 

3 

iv> = 5>>' n '°> 

n=0 

In this case it is easy to declare the object vec directly as: 

Vl:=Matrix(4,2): 
for i from 1 to 4 do 
Vlfrl^A*- 1 : 
Vl[i,2]:=[i-l,i-l,0]: 
end do: 
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However this is not allways the case and if we have a pure state with no known 
f(ni,ri2, ...,71k), but we know which non-zero elements it contains, we can de- 
clare its elements one by one or declare the vec matrix at once. For example to 
declare 

|V>) = |00) + A|11)+A 2 |21) 

we can use: 



d:=2: 




K:=2: 




V:=Matrix(3,2): 


v[i,i] 


=1: 


V[l,2] 


-[0,0]: 


V[2,l] 


=A: 


V[2,2] 


= [1,1]: 


V[3,l] 


=A 2 : 


V[3,2] 


=[2,1]: 


or 



d:=2: 
K:=2: 

V:=Matrix([ fl,f0,0]1, [A,[l,l]], [A 2 ,[2,l]] ]): 



or 



d:=2: 
K:=2: 

V:=«(l, A, A 2 )> | <[0, 0], [1, 1], [2, 1]»: 



Adding StateNormalize (V) for normalization. 

A special family of pure states are readily available in Quantavo: 

Squeezed vacuum: 

for a truncated, unnormalized, pure single mode squeezed vacuum state \<f>) ~ 
Eto^|n)type: 



SqueezedVac(l, d, A); 



for a truncated, unnormalized, pure two mode squeezed vacuum state 
l4>>~£^=oA>,n>type: 



SqueezedVac(2, d, A); 
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These states are given without normalization, since for finite d, \f\ — A 2 doesn't 
normalize them. Calculations and displays are easier this way. 

Coherent States: 

for a truncated single mode coherent state \<p) — J^tZo 77=T \ n ) tyP e: 



Coherent St at e(l, d, a); 



The analytical normalization is rather lengthy so it is left unnormalized. Nor- 
malization can be done at a later time with StateNormalize. 

It is also possible to tensor some vacuum modes to our state. The proce- 
dure TensorVac(-0, s) effectively does the following transformation: \ip) — ► 
\ip) ® |0)® s when applied either to a vec, mat or matcol (see appendix for 
more details). 



Tensor Product: 

The procedure TensorProduct(A, List a, B, Lists) will make the tensor prod- 
uct between modes [List a] of state A and modes [Lists] of state B. So for 
example, 



Co:=CorierentState(l, 4, a); 

Fock:=Matrix([l,[l]]); 

State: =TensorProduct (Co, [1] ,Fock, [2] ) ; 



will result in the following vec object, 

1 [0, 1] 
[1,1] 

1/2 a 2 V2 [2,1] 
_ 1/6 a 3 [3,1] 

Mixed States 



State = 



The generalized states (either pure or mixed) we are interested in can be 
written as: 



d-l 



P = ^2g{n, m)\n)(m\ 



(2) 



where n and fh stand for m, ri2, % and mi, m 2 , ttik respectively. The 
density matrix for these states is d K x d K dimensional which is in general too 
large for the computer to handle. We will therefore describe these states using 
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mat, matcol or poly, objects. 



Two main strategies can be used to declare our states: 

1. If our starting state is pure and will become mixed later, we can declare 
a pure vec object and then convert it into a mat or matcol when needed. 
All procedures to convert are named in the intuitive way: "ob j ect2ob j ect' 
This way to convert vec into mat we have the procedure vec2mat, to 
convert mat into matcol, mat2matcol and so on for all objects and con- 
versions. Therefore, once our pure state vector V has been declared we 
can do the following: 



Vl:=Trim(V): ^eliminate non-zero entries 
Ml:=vec2mat(Vl): ^convert it to a density matrix 



Vl:=Trim(V): 
Mlcol:=vec2matcol(Ml): 



2. if our starting state is mixed to begin with, we can declare our initial state 
as the object mat or matcol. To do so, if we know the functional form 
of g(n,m) in eq.Q then we may directly declare our state. For example 
to declare the state: 

4 

P= Yl A" + ">,n,0)(m,m,0| 

n,m— 

We could use: 



Mlcol:=Matrix(3,25): 
for i from 1 to 5 do 
for j from 1 to 5 do 



Mlcol[i,l]:=A i+J '- 2 ; 

Mlcol[i,2]:=[i-l,i-l,0]; 

Mlcol[i,3]:=[j-l,j-l,0]; 



end do: 
end do: 



Vl:=SqueezedVac(2, 4, lambda): 
Vl:=TensorVac(Vl, 1): 
M 1 col: = vec2matcol ( VI ) : 



or yet again, meaning p - 



E"=oA»)® (£l= A m M|o))] [c.c 



n 



Vl:=SqueezedVac(l, 5, lambda): 

V2:=SqueezedVac(l, 5, lambda): 

V2:=TensorVac(V2, 1): 

V3 : =TensorProduct ( V 1 , [ 1 ] , V2 , [2 , 3] ) ; 

Mlcol:=vec2matcol(V3): 



The basic conclusion is that there is no single way of declaring our state and 
that depending on its structure we have to find a clever way of declaring it. As 
a general guideline, small states without an obvious functional structure can be 
declared giving all elements. Medium sized pure states can be declared using a 
clever loop, 'trimmed' and transformed to vec, mat or matcol. 

It is worth noting that when using the built-in declaration procedures such 
as SqueezedVac, Coherent St ate, TensorVac, Vac, or the inter-converting ones 
like mat2matcol, mat2poly, vec2mat, etc, the values of K and d are automat- 
ically updated. However, if the states are declared from scratch, the values 
of K and d should be explicitly declared. Applying to our state p the proce- 
dure f indKnd(p) can help troubleshoot by reevaluating the value of these global 
variables. 

3.2.2 Evolution 
Beam Splitter 

We can apply a beam splitter to a vec or matcol state V with the BS procedure: 



Vl:= BS(V,i,j): 



Where i, j specify which modes the beam splitter acts on. It is therefore 
essential to carefully label the modes of our quantum optical circuit. This will 
effectively do the following mode transformation: 

Leaving T and V as unevaluated variables. If we want to use different re- 
flectivities and transmittivities for our Beam Splitter (variables or numbers) we 
may use myBS and input H = to' and 'r = r ' as follows: 



Vl:= myBS( V , i , j , t , r ): 



For the mixed state object matcol we may use the same procedures and the 
effective transformation will be: 
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Further details about its use can be found in the example in section {I} 
Phase Shifter 

The phase shifter procedure (PS) can be used with mixed and pure states: 



Vl:= PS( vec/matcol , i ,0 ): 



The input vec/matcol means that either vec or matcol objects can be given 
as inputs. This procedure makes the effective transformation: 

a[ = e^cti (6) 
K = (7) 



Build Unitary 

If we wish to construct a unitary matrix that transforms the modes of light 
and describes a given linear optics (LO) circuit we may use the procedure 
BuildUnitary. Together with UnitaryEvolutionit will evolve the state through 
a given LO circuit. To build a unitary consisting of Beam Splitters (BS) and 
Phase Shifters (PS) we will do the following: 



U:= BuildUnitaryQList]): 



This will create a matrix U of dimension K x K, that later will transform 
the K modes. The list that is BuildUnitary's input must have a precise for- 
mat. It must be a list of lists. For example: List:=[[l,2,t,r],[3,phi],[3,4,q]]; 
means that first a BS with transmitivity t 2 and reflectivity r 2 will be applied 
to modes 1 and 2, then a PS will be applied to mode 3 and finally a BS with 
transmitivity q 2 and reflectivity (1 — q 2 ) will be the last operation. If we have 
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4 modes, this will build the matrix 



U = 



( t r 

-r t 



\ 

/ 1 o 
1 



V o o 








o 





Q 



/ 



\ / 1 \ 

10 



\ o o o i j 



( t r \ 

-r t 

10 

\ 1 j 



In general, in our list of lists, lists with 4 elements, like [i, j, t, r] build BS 
transformations between modes i and j, lists with 3 elements like [i, j, t] build 
BS transformations for modes i and j such that r = vl — t 2 , and lists with two 
elements like [i, <p] build PS transformations on mode i. 

UnitaryEvolution 

Whether we have just built a unitary matrix with BuildUnitary or we have 
a K x K arbitrary unitary matrix to transform our modes {ai}, we can use this 
procedure as follows: 



V:= UnitaryEvolution(U, vec/matcol): 



where U is the unitary matrix of dimension K x K and our state is described 
by a vec or a matcol object. This will effectively implement the mode trans- 
formation: 

a! = Ua 

3.2.3 Measurement 

The Project Procedure: 

If we wish to know the state after a projective measurement we may use the 
procedure Project. Depending on the inputs we give to the procedure it will do 
a projective measurement and return a density matrix (matcol) or state vector 
(vec). Below is a description of its different uses: 

1. if given (veci, list, vec 2 ) and say veci and vec 2 describe respectively 
\ipi) and \ip2), Project(veci, list, vec 2 ) returns the vec (in principle 
unnormalized) corresponding to the expresion: 



l^2> = {\lpl)(lpl\list ® Irest) |^2> 
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As an example consider list=[2,3] meaning that we want to measure modes 
two an three. The vec object that corresponds to the projector 
must therefore have kets with 2 modes. 

For example, if ^2) is 



and |V>i) is 



Then, 





1 


|0000 > 


V2 := 


X 


|1100 > 






|2200 > 


VI := 


= [1 


1 10 > ] 



S := Project(Vl, [2,3], V2); 



will return: 



[ 



11100 > 



2. if given (matcol 1, list, vec 2) and say matcol 1 and vec 2 describe M n 
and ^2) then it will return the matcol object: 

P = ( M nn st ® Ire-t) \^h)(H (M„ Hst ® l rest ) f 

So for example, taking the same V2 as above and the Kraus Operator (or 
projector): 



M := 



1 [1,1] [1,1] 
1 [0,0] [0,0] 



Then, 



S := Project(M, [1,2], V2); 



will return the state: 



S = 



which is the density matrix corresponding to: 

1 [0,0,0,0] 



1 


[0,0,0,0] 


[0,0,0,0] 


X 


[0,0,0,0] 


[1,1,0,0] 


X 


[1,1,0,0] 


[0,0,0,0] 


XX 


[1,1,0,0] 


[1,1,0,0] 



x [1,1,0,0] 
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3. if given (vec 1, list, matcol 2) and say vec 1, matcol 2 describe \ipi) and 
p respectively, then it will build IV'iXV'il and return the matcol object 
corresponding to: 

P' = {\^l){^l\list ® Irest) P (\lpl)(^l\ust ® lrest) 



4. if given (matcol 1, list, matcol 2) and say matcol 1 and matcol 2 

describe M n and p respectively, then it will return: 

p' = {M nust ® l rest ) p (Mn,,,, <g> I rest ) f 

(see further down for POVM measurements) 
In a nutshell: 

Project(|V>i), list, |t/> 2 » -> (IV'lXV'llHst ® lrest) 1^2) 

Project(M„, list, IV2)) -» (M„ H , t ® I rest ) I^KV^I (K lj8l ® lrest) 1 

ProjeCt(|^i), list, p ) -> (IV'lKV'llzist ® lrest) P (IV'l) (^1 1 Jist ® lrest) 

Project(M„, list, p ) -» (M„ H , t ® I rest ) p (M„ !ist ® l^) 1 

If we encounter a destructive measurement, it is possible to trace out the mea- 
sured modes with the Traceout procedure (cf. appendix). Also POVMresult 
traces out the measured modes (see next section). We may otherwise calculate 
the full trace of a mat or matcol object with StateTrace or multiply vec and 
matcol in different orders thanks to StateMultiply. 

The Probability Procedure: 

This procedure will calculate the probability of a measurement result or an 
expected value. It considers the same cases and objects as the above Project 
procedure. It uses the definition, 

Tr {E n p} 
Tr{p} 

and assumes that {E n } satisfy E n > and ^2 n E n = II to calculate probabil- 
ities. Therefore, one should verify that these conditions hold for the matrices 
describing E n in order to obtain meaningful probabilities. Below we show more 
details for different inputs: 

Pro^imydW, list, | W ) MlMj^i^ KWM 

Tr {|r/> 2 ) (ip 2 1} 

where (IV'l), IV^}) are converted to matcol objects in an intermediate step. 
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One should pay attention to the choice of the projection operator |^>i)(V>i|. If 
it is not nomalized it can give unphysical values for the probability. 

Now for a given POVM or Projector, 

Probability^, hst, |V> 2 » - 

Or given a vec and a density operator matcol: 

Probability^), hst, p ) -» — -— , — 

Tr {p} 

where IV'i) has been converted to a matcol object. 



And finally for a POVM and a density operator: 

Tr {(E n list ® Irest) P } 



Probability(i?„, list, p ) 



Tr{p} 



3.2.4 POVM measurements 

Quantavo posesses a procedure to describe POVM measurements. If our state 
before the measurement is \ip) or p, and the POVM elements are described by 
the set {Em}, satifying ^ m E m — II and E m > 0. Then the state after the 
measurement will be: 

p' = Tr iij; .. fc {E m p} 

or 

p / = Tr iij; .. fc {£ m |V>)W>|} 
assuming Tr {p} = 1 or Tr {|V')('^|} = 1- 

This will be implemented by the procedure POVMresult which will take as in- 
puts, 

POVMre suit (matcol, List, matcol/vec) 

and implement the operation, 

POVMresult(£ m , List, p) = 

Tr{p} 

where 'rest' are all the indexes not included in l List\ 



3.2.5 Declaring POVMs 

The matcol objects that represent POVM operators need to be declared. One 
option is to declare them as standard states and convert them to matcol. Since 
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avalanche photo diode detectors (APDs) are a standard tool in quantum optics, 
an interactive tool to declare them is also provided. Executing, 



APDQ; 



Will bring up an interactive menu, where we can choose if we want the |0)(0| or 
I — |0)(0| event, how many photons we will consider and what loss should be 
added to it. For example choosing the input (0, r, 4) will return, 



1 


[0] 


[0] " 


r 2 


[1] 


[1] 


r 4 


[2] 


[2] 


r 6 


[3] 


[3] 


r 8 


[4] 


[4] . 



[1] " 

[2] 
[3] 
[4] . 

recovering the expected tto + tti = I. Note the convention for the BS in front of 
the detector for which r = is a perfect detector and r 2 + t 2 = 1. 

These are therefore the main tools to describe and simulate quantum mea- 
surements in this framework. 

3.2.6 State Properties 

Some questions that will interest us will concern hermiticity, normalization and 
entanglement measures. So far, Quantavo has the following usefull procedures: 

IsHermitian to check for selfadjointness, 
IsNormalized to check for normalization, 
StateNormalize to normalize mat, matcol or vec 
StatePartialTranspose to partial transpose p T 

Negativity to calculate the negativity as the sum of all negative eigenvalues of 

the partial transposed density matrix. 

LogNegativity to calculate the Logarithmic Negativity. 

StateApprox to do symbolic or numeric approximations transforming our mat- 
col or vec state. 

For practical examples on how to use them one can refer to the following 
section. Otherwise, a detailed dictionary of procedures can be found in the 
appendix. 



and choosing (1, r, 4) will return, 



" 1-r 2 


[1] 


1 -r 4 


[2] 


1-r 6 


[3] 


1-r 8 


[4] 
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4 Practical Example 

4.1 Squeezed state photon subtraction 




Figure 1: Setup of the photon subtraction. 

Our initial state is a pure two mode squeezed state that can be described 
by IV'a) = Vl — A 2 YlnLo ^™l n : n )- We would like to declare it and propagate 
it through the circuit presented in fig. Q] We observe that our initial state has 
three modes, one of which is a vacuum mode. We can work out as an example 
the state with up to 4 photons: 



# Create a truncated two mode squeezed vacuum state 

V := SqueezedVac(2, 5, lambda); 

# Then add the vacuum mode: 

V := TensorVac(V,l); 



This will output state: 



1 


[0,0,0 


A 


[1,1,0 


A 2 


[2,2,0 


A 3 


[3,3,0 


A 4 


[4,4,0 



We then apply the corresponding Beamsplitter transformation (if "d" changes 
from the BS transformation, it is automatically recalculated after the BS oper- 
ation and reset to its new value. In this case it doesn't change). 



V1:=BS(V,1,3); 



which returns, if displayed with "Display State", that is Dstate(Vl), the fol- 
lowing output: 
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VI 



1 

xt 

X 2 t 2 

A 3 t 3 

A 4 t 4 

Xr 

V2X 2 tr 

V3A 3 t 2 r 

2 XH 3 r 
X 2 r 2 

VSXHr 2 

V6XW 
A 3 r 3 

2 A 4 tr 3 
A 4 r 4 



Measurements: 



'| ooo >'■ 
'jllO >'■ 
'|220 >'■ 
'|330 >'■ 
'|440 >'' 

ion >'■ 

'| 121 >'' 
'|231 >'' 
'|341 >'■ 
'|022 >'■ 
'|132 >'■ 
'|242 >'■ 
'|033 >'' 
'| 143 > h 
'|044 >'■ 



First we will consider the case of a perfect photon number resolving detec- 
tor. The projector |1)(1| can simply be introduced as: 



Proj:=Matrix([l,[l]]); 



And the state after the measurement will be obtained with: 



V2:=Project(Proj,[3],Vl); 



V2 = 



returning: 

Xr [0,1,1] 

V2X 2 tr [1,2,1] 

V3X 3 t 2 r [2,3,1] 

2A 4 i 3 r [3,4,1] 

We may now trace-out the measured mode. Tracing out can then be done as 
follows: 
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Ml:=Traceout(V2,3); 



We may also simulate a measurement with an avalanche photo-diode de- 
scribed by the Kraus Operator: 6 = 1- |0)(0|. That is, if the state "VI" 
represents the state p, we want to find the state p' resulting from the measure- 
ment: 

ft ~ (o 3 ® ii, 2 )/o(6 3 ® ii, 2 ) t 

For this simple example we have at most 4 photons so we will approximate 
6 ~ |1)(1| + |2)(2| + |3> <3 1 + |4)(4|. Therefore to measure mode 3 wc construct 
the associated POVM which can then be expressed as a matcol object: 



POVM := Matrix(4, 3); 

for i from 1 to 4 do 

POVM[i,l]:=l: 

POVM[i,2]:=[i]: 

POVM[i,3]:=[ij: 

od: 



And the state after a 'click' in the detector (tracing out this inaccessible mode) 
will be: 



M2:=Project(POVM,[3],Vl): 
M3:=Traceout(M2,3): 



or directly 



M3:=POVMresult(POVM,[3],Vl): 



Approximations: 

M3 has 30 coefficients. We may want to use the approximation procedures 
to simplify our calculations. If for example wc know that r << 1 and A << 1 
we may wish to delete all terms in the density matrix for which the coefficients 
containing r "A m satisfy n+m > 7 (therefore getting rid of all the small terms up 
to the chosen order). The procedure StateApprox can be used for this purpose. 
This way, 
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M4: 


=RealDomain:-simplify (M3) : 


# assume real variables for simplicity 


M5: 


=StateApprox(M4, [lambda, r 


,7); 



Will deliver 



M5 = 



A 2 r 2 


[0,1] 


[0,1] " 


V2\Hr 2 


[0,1] 


[1,2] 


\H 2 r 2 V3 


[0,1] 


[2,3] 


2 X 5 t 3 r 2 


[0,1] 


[3,4] 


V2\Hr 2 


[1,2] 


[0,1] 


2 \H 2 r 2 


[1,2] 


[1,2] 


V6X 5 t 3 r 2 


[1,2] 


[2,3] 


\H 2 r 2 V3 


[2,3] 


[0,1] 


\/6AW 


[2,3] 


[1,2] 


2 AW 


[3,4] 


[0,1] . 



See the dictionary in the appendix for details on how to use the StateApprox 
procedure. 

Entanglement 



We may now compute different properties of this state as for example the Neg- 
ativity. We take the state after the detection of exactly one photon and obtain: 



Neg:=simplify (Negativity (Ml)) assuming 


t::positive, 


r::positive, 


lambda: :positive; 







which will return: 



At (2 \ 2 t 2 + V2 + \/3A t + V3X 2 t 2 V2 + 2 \H 3 V2 + 2 \H 4 V3) 
NeQ ~ 4A 6 t 6 + 3A 4 i 4 + 2X 2 t 2 + 1 

The eigenvalues can be hard to find if our state has symbolic entries with com- 
plex conjugation. A previous simplification (assuming real variables) can help, 
but other solutions are possible. This explicit expression allows us to plot dif- 
ferent values of the negativity for various ranges of parameters: 



plot3d(Neg(A,t,\/l - < 2 ), A=0..1,t=0..1, axcs=boxed); 

will display the plot in fig O: 

Further details, examples and documentation can be found on the website 
of the Imperial College London quantum information group, 
http : //www . imperial . ac .uk/quantuminf ormat ion/research/downloads 
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Figure 2: Logarithmic Negativity (t, A). 
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A Dictionary of procedures: 

A quick reference list is provided and an alphabetically ordered dictionary fol- 
lows. 



Objects 



Appearance 



Declaration 



vec 



1 


[0,0,0] 


A 


[1,1,0] 


A 2 


[2,2,0] 


A 3 


[3,3,0] 



vcc:=CoherentState(K,d,a) 

vec:=SqueezedVac(K,d,A) 

vac:=Vac(K); 

vecl:=TensorVac(vcc,m) 

vecl:=TensorProduct(V,[l,2],W,[3,4]) 

vec:=Trim(V) 



mat 



[0,0,0] [1,1,0] 
[0,0,0] 1 A 
[1,1,0] A A A 



mat:=vec2mat(vec) 
mat:=matcol2mat (matcol) 
or direct declaration of the 
Matrix 



matcol 



Pmatcol 



1 


[0,0,0] 


[0,0,0] 


A 


[0,0,0] 


[1,1,0] 


A 


[1,1,0] 


[0,0,0] 


AA 


[1,1,0] 


[1,1,0] 


A 2 A 2 


[2,2,0] 


[2,2,0] 



matcol:=vec2matcol(vec) 
matcol: =mat 2mat c o 1 (mat ) 
or direct declaration of the 
3 column Matrix 



CoherentState(if , d 7 a) 

ready made states are: SqueezedVac(i<:, d, A) 
TensorVac(vec, to) 

IdentityState(d,K) 



d-1 

£ 

n=0 



n) 



£a» 



n=0 



10) - \4>) ® w 



«-d K xd" 
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Common Procedures: 


Linear Optics: 


Measurements 


BS(vec/matcol, i, j) 




myBS(vec/matcol, i, j, t, r) 


Project(vec/matcol, [List], vec/matcol) 




Probability( vec/matcol, [List] vec/matcol) 


PS(vec/matcol, i, <j>) 


Traceout(mat/matcol, i) 




POVMresult(matcol, [List], vec/matcol) 


BuildUnitaryQLzsi of Lists]) 


APD({0,l},r,d+l) 


UnitaryEvolution([7, vec/matcol) 





Display 


Algebraic Operations 


Dstate(vec/matcol/mat) 
PlotState(vec/matcol/mat, width, height) 


StateApprox(vec/matcol, list, N) 

St at eMult iply(vec / matcol, vec / matcol) 

StateComplexConj ugate(mat/matcol) 

StateNorm(vec) 

StateNormalize(vec/mat / matcol) 
StateTrace(mat /matcol) 
Traceout(mat / matcol, i) 

TensorProduct(vec/matcol, list, vec/matcol, list) 
StatePartialTranspose(mat/matcol, i) 
deltaK(iJ) 




State Properties: 


Entanglement & Energy 


IsHermitian(matcol/mat) 
IsNormalized(vec / matcol/mat) 
StateNorm(vec) 
StateSort(vec/matcol) 
FindKnd(vec / matcol/mat) 


Negat ivity(vec / matcol/mat) 
LogNegativity(vecmatcol/mat) 
Entropy(vec/matcol) 
Energy(vec /matcol) 
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A 



APD(0/l,r,N): This procedure is interactive and is called by executing, 
APDQ; 

It declares the POVM describing a lossy avalanche photo diode detector (APD). The 
input "0" or "1" will select between the no-click or click events respectively, "r" 
will be the amplitude r of the reflectivity R = r 2 of the BS in front of the detector 
characterizing its loss, and N will be the maximum number of photons. For example 
choosing the input (0, r, 4) will return, 



" 1 


[0] 


[0] " 


r 2 


[1] 


[1] 


r 4 


[2] 


[2] 


r 6 


[3] 


[3] 


. r 8 


[4] 


[4] . 



and choosing (1, r, 4) will return, 



- ! _ r 2 


[1] 


[1] " 


1-r 4 


[2] 


[2] 


1-r 6 


[3] 


[3] 


.1-r 8 


[4] 


[4] . 



recovering the expected 7To + 7Tl = J. 



[ Input: ({0,1}, r £ [0,1], positive integer), Output: Matrix, Calls: Quantavo, Lin- 
ear Algebra] 

B 

BS(vec/matcol, i, j): This will effectuate the Beam Splitter transformation: 




for the chosen modes. That is eq. (8) for objects of type vec and eq. (8) and (9) for 
objects of type matcol. The value of '<f and 'K' are evaluated and reset to the actual 
value after this operation. 

[ Input: (Matrix, integer, integer), Output: Matrix, Calls: Quantavo, LinearAlgebra, 
poly2matcol, matcol2poly, matcolBS, vecBS] 

BuildUnitary(list of lists) : This procedure will build a K x K unitary matrix 
corresponding to a linear optics circuit consisting of beam splitters (BS) and phase 
shifters (PS). To do so a list of lists must be provided. The lists inside the list can 
have 2, 3, or 4 elements, and should be in the same order as we want to apply those 
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transformations in the circuit. Lists with 2 elements will be considered PS and lists 
with 3 or 4 elements as BS: 

The list [i, £] will implement a PS on mode i. 

The list [i,j,t] will implement a BS on modes i and j with transmittivity t 2 = T and 
reflectivity R = 1 - T = 1 - t 2 . 

The list [i,j,t,r] will implement a BS on modes i and j with transmittivity t 2 = T and 
reflectivity r 2 = R. 

To implement the transformations one after another we can give for example the 
input 

list = [[i u £],[i 2 , 32, t],[i 3 , h, t, r],[i 3 ,0]]. 

[ Input: List, Output: Matrix, Calls: Quantavo, LinearAlgebra] 

c 

CoherentState(m,d,a): This procedure builds an object of type vec that de- 
scribes the state \<f>) = X) n =o l n )' 8m - It is not normalized 

[ Input: (whole number, whole number), Output: Matrix, Calls: Quantavo, LinearAl- 
gebra] 

D 

DP(Matrix, Matrix) : This procedure makes the Kronecker/Direct /Tensor prod- 
uct between any two matrices. As long as it is declared as a matrix, it can also make 
the Kronecker product between vectors and matrices. The vector, however would have 
to be declared as a Matrix. For example as V := Matrix([[l], [2], [0], [a]]). 
[ Input: Matrix, Output: Matrix, Calls: Quantavo, LinearAlgebra] 

DeltaK(i, j): This procedure takes the Kronecker Delta from two inputs. It works 

both symbolically and numerically. For example: 

deltaK(l,l)=l, deltaK(l,2)=0, 

deltaK(a,a) = l, (even if a is not specified). 

deltaK(a,b); will return Quantavo:-deltaK(a, b) 

However, if we execute "a:=b;" and evaluate it again it will be 1 

[ Input: (string or number, string or number), Output: 0,1, unevaluated string, Calls: 
Quantavo] 

Dbra(List): This is a Display procedure. When given a list which stands for a 
bra, for example List =[0,1,3] it builds a bra for display and returns (013 1 . 
[ Input: List, Output: ( string|, Calls: Quantavo]. 

Dbraket(List, List): This is a Display procedure. When given two lists which 
stand for a ket and a bra it returns a \ket){bra\. for example IList=[0,l,l] and JList 
= [0,1,3] then Dbraket(IList, JList) =|011)(013|. 
[ Input: (List, List), Output: | string) ( string|, Calls: Quantavo] 

Dket(List): This is a Display procedure: When given a list which stands for a 
ket, for example List =[0,1,3] it builds a ket for display and returns |013). 
[ Input: List, Output: J string), Calls: Quantavo] 

Dstate(vec/mat/matcol): This is a Display procedure. When given ei- 
ther a mat, vec or matcol object it transforms the lists of modes into bras and kets. 
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For example if it is given the object mat: 

[0,0,0,0] [1,1,0,0] 

p= [0,0,0,0] 1 A 

[1,1,0,0] A A A 
it will display bras and kets in the following way: 

< 0000| < 1100 

Dstate(p) 



1 0000 > 1 
11100 > A 



A 
AA 



[ Input: Matrix, Output: Matrix, Calls: Quantavo, LinearAlgebra] 

E 

Energy(vec/matcol): This will output the expected value of the energy of 
our state defined as, 

(E) = {i>\E\i>) 

with 

(E) = ?^Tr{(a 1 a + l/2)p} 

or for multipartite states, 

(E) = huTr {(iVi + iV 2 + ... + N K + K/2)p} 

[ Input: Matrix, Output: Number or Analytic Expression, Calls: Quantavo, LinearAl- 
gebra] 

Entropy (vec/matcol): This will output the entropy of our state defined as, 

S = ^2 A l /og 2 (A l ) 

where A; are the eigenvalues. The procedure transforms the objects vec or matcol 
into a mat density matrix and then finds the eigenvalues. 

[ Input: Matrix, Output: Number or Analytic Expression, Calls: Quantavo, LinearAl- 
gebra] 

EvalState(vec/matcol): This simply applies the Maple function evalf to 
the first column of our vec or matcol objects. 

[ Input: Matrix, Output: Number or Analytic Expression, Calls: Quantavo, LinearAl- 
gebra] 



f indKnd(vec/matcol/mat): This will search through the given state to find 
the number of modes K and the dimension of the modes d. It updates these global 
variables with the K and d found and displays them as output 
[ Input: Matrix, Output: (K,d), Calls: Quantavo, LinearAlgebra] findKnd 



I 
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IdentityState(Nr of photons, Nr of modes): This will declare an 
identity matrix represented as a matcol object. For example, for (2,2) as input we 
will obtain 



1 


[0,0] 


[0,0] ■ 


1 


[0,1] 


[0,1] 


1 


[0,2] 


[0,2] 


1 


[1,0] 


[1,0] 


1 


[1,1] 


[1,1] 


1 


[1,2] 


[1,2] 


1 


[2,0] 


[2,0] 


1 


[2,1] 


[2,1] 


1 


[2,2] 


[2,2] . 



[ Input: (integer, integer), Output: Matrix, Calls: Quantavo, LinearAlgebra] 

indexstate (vec / matcol) : This procedure will transform the last two columns 
of a vec or matcol object. Each list containing the number of photons in each mode 
will be translated into a natural number with the procedure VectorRow(List,d). This 
number will indicate the order in which the modes are ordered. In effect it does a 
number basis change from a d-base to a 10-base. For example, 



" 1 


[o,o] - 




I I 




[1,1] 




£ 5 




[2,2] . 




_ e 9 _ 



[ Input: 3 column Matrix, Output: 3 column Matrix, Calls: Quantavo, LinearAlgebra] 

IsHermiti an(matcol/mat): This procedure verifies if its input describes 
a Hermitian state. Therefore if the density matrix it represents verifies p + — p. If 
it is Hermitian the value returned will be 0, otherwise it will be 1. Maple might not 
recognize products of conjugated complex variables as equal, so one has to make sure 
they are simplified. 

[ Input: Matrix, Output: or 1 and printed answer, Calls: Quantavo, LinearAlgebra] 

IsNormalized(vec/matcol/mat): This procedure checks if the considered 
object is normalized. Therefore for vec if (ip\ip) = 1, for matcol and mat if Tr(p)=l. 
If the state is not normalized it will return its norm or trace. 
[ Input: Matrix, Output: printed answer, Calls: Quantavo, LinearAlgebra] 

L 

LogNegat ivity(vec / matcol / mat) : This procedure evaluates the Logarith- 
mic negativity according to: 

LogNegativity(M) := log 2 (2 * Negativity(M) + 1) 

See Negativity for further details. 

[ Input: Matrix, Output: Expression or number, Calls: Quantavo, Negativity] 

M 
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matcol2ma~t(matcol): This translation procedure takes an object of type mat- 
col and converts it into an object of type mat. To do so, it sorts the matcol object, 
counts the number of different bras and kets needed and constructs the hermitian ma- 
trix mat with the least number of zeros containing all the elements of matcol 
[ Input: Matrix, Output: Matrix, Calls: Quantavo, LinearAlgebra, StateSort] 

matcol2poly(matcol) : This translation procedure takes an object of type mat- 
col and converts it into an object of type poly. To do so it looks at the |ket) and (bra| 
of each row of the Matrix and builds the associated monomial with its corresponding 
coefficient. 

[ Input: Matrix, Output: Polynomial, Calls: Quantavo, LinearAlgebra] 

mat2matCOl(mat) : This translation procedure takes an object of type mat and 
converts it into an object of type matcol 

[ Input: Matrix, Output: Matrix, Calls: Quantavo, LinearAlgebra, Trim] 

Diat2poly(mat) : This translation procedure takes an object of type mat and 
converts it into an object of type poly 

[ Input: Matrix, Output: polynomial, Calls: Quantavo, LinearAlgebra] 

modesmatCOl(matcol with numbers): This procedure is the inverse of 
the procedure indexmatcol. It will transform a 3 column matrix, that has numbers in 
the 2nd and 3rd column to one that has the equivalent modes (lists with the number 
of photons) 

[ Input: 3 column Matrix, Output: 3 column Matrix, Calls: Quantavo, LinearAlgebra] 

myBS(matcol/veC, i, j, t, r): This procedure is the same as BS, but it 
gives the user the option to choose the transmittivity T — t 2 and reflectivity R = r 2 
of the beam splitter. 

[ Input: Matrix, Output: Matrix, Calls: Quantavo, LinearAlgebra, poly2matcol, mat- 
col2poly, mymatcolBS, myvecBS] 

N 

Negativity(vec/matcol/mat): This procedure finds the eigenvalues of the 
partial transposed density matrix, and calculates the sum of all the corresponding 
non-negative eigenvalues (it also divides by the Trace of the state in case it is not 
normalized): 

Negativity = ^ N_A 

where \i are the eigenvalues of the partially transposed density matrix 

[ Input: Matrix, Output: Expression or number, Calls: Quantavo, LinearAlgebra, 

StatePartialTranspose] 

P 

PlotState(vec/matcol/mat, W, h): The coefficients of our states need 
to be real numbers before we can plot the state vector or density operator. 

The procedure will plot a vec object as a bar diagram. In the abscissa, all possi- 
ble kets between |0, 0, .., 0} and \(d — 1), (d — 1), (d — 1)) are labeled from 1 to d K . 
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For example, 



1.0 


|0000 > " 


0.5 


1100 > 


0.25 


|2200 > _ 



with K=2 and d=3, will be displayed as shown in Fig. [3] 




2 4 6 s 10 



Figure 3: Plot for a pure state (vec object). 

To plot a matcol object, the width w and height h of the bars in the bar diagram 
must be given as inputs (w=0.5 is the maximum width for the columns not to overlap). 
A 3D plot will be the output. Bras and kets will also be labeled from 1 to d K . The 
density matrix for the above state (w=0.5, h=10) is shown in Fig. U] 




Figure 4: Plot for a density matrix (matcol/mat object). 

[ Input: (Matrix, width, height), Output: , Calls: plots, geom3d, Quantavo, LinearAl- 
gebra] 

poly2mat COl (poly): This procedure transforms a polynomial of the modes 

d d K 

Poly = /(n, in) J^J a" 3 b™ 1 into an object of type matcol 

mi ,1712 !■■ n li Tt 2 >--i n K J— 1 
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[ Input: polynomial, Output: Matrix, Calls: Quantavo, LinearAlgebra] 



poly2vec(poly): This procedure transforms a polynomial of the modes 

d K 

Poly = f(n) Y\ a j J ^ nto an °bject of type vec 

n 1 ,n 2 ,..,n K j — 1 

[ Input: polynomial, Output: Matrix, Calls: Quantavo, LinearAlgebra] 

P0VMreSUlt(matCOl, List, vec/matcol): If our state before the mea- 
surement is \ip) or p, and the POVM elements are described by the set {E m } then the 
unnormalized state after the measurement will be: 
POVMresult(_B m , List, p) = Tr Llst {(E mList <g> lrest) p} 
where l rest' are all the indexes not included in 'List'' '. 

[ Input: Matrix, List, Matrix, Output: Matrix, Calls: Quantavo, LinearAlgebra] 

Probability(vec/matcol, List, vec/matcol): This procedure will cal- 
culate the probability of a measurement result or an expected value. It considers the 
same cases and objects as the Project procedure. It uses the definition, 

Tr {E n p} 
Tr{p} 

and assumes E n > and J] n E n = 1 to calculate probabilities. Therefore, one should 
verify that these conditions hold for the matrices describing E n in order to obtain 
meaningful probabilities. Below we show more details for different inputs. 

Probability Vi , list, h/>2 -> — rr - 

Tr {|V>2>(^2|} 

where (\ipi), \ip2)) are converted to matcol objects in an intermediate step. One 
should pay attention to the choice of the projection operator \ipi){if)i\. If it is not 
nomalized it can give unphysical values for the probability. 

Now for a given POVM or Projector, 

where E n is assumed to be a matcol object. 
Or given a vec and a density operator matcol: 

Probability( Vi), list, p ) -> — — ' f , — 

where \tpi) has been converted to a matcol object. 



And finally for a POVM and a density operator: 

Tr {(E n ® I rest ) p } 



Probability (E„, list, p ) 



Tr{p} 



[ Input: (Matrix, list, Matrix), Output: Expression or number, Calls: Quantavo, Lin- 
earAlgebra, StateTrace, StateNormalize, vec2matcol] 

Project (vec/matcol, list, vec/matcol): This procedure returns the 



32 



state after a measurement in the following cases: 



1. if given (veci, list, vec2) and say veci and vec2 describe respectively 
and \ip2), Project(veci, list, vec2) returns the vec (in principle unnormalized) 
corresponding to the expresion: 



W2) = (|-0i><^i|ii»t ® lr 



Note that if for example list=[l,2,3] it means that we want to measure modes 
one, two an three. The vec object that corresponds to \ipi) must therefore have 
kets with 3 modes. 



For example, if \ip2) is 



and \ipi) is 



V2 := 



x 

~2 



VI := 



1 0000 > 
1 1100 > 
1 2200 > 

|10 > ] 



Then, 



S := Project(Vl, [2,3], V2); 



will return: 



S=[ x 1100 > ] 



2. if given (matcol 1, list, vec 2) and say mated 1 and vec 2 describe M n and 
|^2) then it will return the matcol object: 

P = ( M ™list ® ^rest) I^X^I (M„ Kst (g) IrestY 



So for example, taking the same V2 as above and the POVM: 

1 [0,0] [1,1] 

M := 



1 [1, 1] [0, 0] 
1 [1,1] [1,1] 



Then, 



S := Project(M, [1,2], V2); 
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will return the state: 

xx [0,0,0,0] [0,0,0,0] 

x + xx [0,0,0,0] [1,1,0,0] 

xx + x [1,1,0,0] [0,0,0,0] 

x + xx + x + 1 [1,1,0,0] [1,1,0,0] 
which is the mixed state that corresponds to: 

x [0,0,0,0] 



5' 



x + 1 [1,1,0,0] 



3. if given (vec 1, list, matcol 2) and say vec 1, matcol 2 describe \ipi) and 
p respectively, then it will build \ipi) (ipi | and return the matcol object corre- 
sponding to: 

P = (|^l)(V , l[«st 8> Irest) p (\lpl){lpl\list ® Irest) 



4. if given (matcol 1, list, matcol 2) and say matcol 1 and matcol 2 describe 
M n and p respectively, then it will return: 

p = (M„ ljs{ ® Irest) p (M„ list ® l r est) + 

In a nutshell: 

Project(|V>i), list, \ip2}) 
Project(Af n , list, \tp 2 )) 
Project(|?/>i), list, p ) 
Project (M n , list, p ) 

[ Input: (Matrix, list Matrix), Output: Matrix, Calls: Quantavo, StateComplexConju- 
gate, LinearAlgebra, vec2matcol] 

PS(vec/matcol, j, (j)): PHASE SHIFTER: This procedure makes a phase shifter 
transformation to our state. It is implemented making the following mode transfor- 
mation to the specified mode j: 

a'j — e^aj (10) 
b'j = e-'% (11) 

(that is transformation (|10|) for a vec and transformation (|10|) and (|lip for matcol 
and mat.) 

[ Input: (Matrix, whole number, symbol or number), Output: Matrix, Calls: Quantavo, 
LinearAlgebra] 

s 

SqueezedVac(m,d,A): This procedure builds a vec describing the state \<fr) ~ 
X^n=o ^™\ n )® m where m can be m = 1, 2. The state is not normalized. 



— * (\lpl) (tpl\list ® Irest) \i>2) 

-> (M„ Ust ® Irest) I^X^I (M„ |js , ® Irest) 1 " 

-> (|V>l) (^1 | fist ® Irest) |0 ( |V>l) (V>1 |«st ® Irest) 

-» (M„„ jf g> l re st) P (M„ Kjl ® lrest) + 
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[ Input: (whole number, whole number, string or number), Output: 2 column Matrix, 
Calls: Quantavo, LinearAlgebra] 



StateApprox(vec/matcol, list, n): This procedure is used to reduce the 
size and complexity of the vec or matcol objects with an approximation. 

Numerical Approximation: If the coefficients of our state are numbers then use as 
follows: StateApprox(vec/matcol, [], n). All Rows for which the coefficient < I0~™ 
will be deleted. 



Symbolic Approximation: If the coefficients are symbolic polynomials (or monomi- 
als) and some variables are small, we may choose to delete all terms with a certain 
power in those variables. If we execute S := StateApprox(M,[Zam6do,r],10), then all 
terms containing \ n r m such that n + m > 10 will be deleted. For example if we have 
a state of the form: 



M 



I + y 4 x 2 + yx 4 
y 7 x + x 3 + x 5 
x 2 + x 4 + x 6 + y 5 



[0,0,0,0] 
[1,1,0,0] 
[2,2,0,0] 



then, 



S := StateApprox(M,[y,x],5); 



will return the state: 



1 + yx 4 



[0,0,0,0] 



x a + x 5 [1,1,0,0] 
x 2 + x 4 + y 5 [2,2,0,0] . 



and 

S := StateApprox(M,[y,x],f); 



will return the state: 

[ 1 [0,0,0,0] ] 

[ Input: (Matrix, list, integer), Output: Matrix, Calls: Quantavo, LinearAlgebra] 

StateComplexConjugate(matcol/mat): It will complex conjugate all 

the coefficients of an object of type mat or matcol. 

[ Input: Matrix, Output: Matrix, Calls: Quantavo, LinearAlgebra] 

StateMultiply(number/matCOl, matcol/vec): This procedure mul- 
tiplies following the matrix multiplication rules the following objects: 

k x matcol, 
k x vec, 

matcol x matcol 
matcol x vec 
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where k is a number, a variable, a function or a string 

[ Input: (string/Matrix, Matrix), Output: Matrix, Calls: Quantavo, LinearAlgebra] 

StateSort(vec/matcol): This procedure sorts by rows the objects vec or 
matcol. The order is determined in the following way. The modes are converted from 
the numeral base d to the numeral base 10 and are then sorted by increasing order. 
For matcol they are sorted first according to the order of the "kets" and inside the 
same "ket" number, by the order of the "bras". For example, with d — 2 and K = 2, 
[A, [0, 0]] < [A 2 , [1, 0]] < [A, [0, 1]] < [A 3 , [1, 1]]. The output is the ordered state. 
[ Input: Matrix, Output: Matrix, Calls: Quantavo, LinearAlgebra, Tribullesmatcol, 
Tribullesvec] 

St at eTrace (mat /matcol): This procedure performs the trace on mat or 
matcol states. It simply adds all diagonal elements. 

[ Input: Matrix, Output: Real Number or symbolic expression, Calls: Quantavo, Lin- 
earAlgebra] 

StateNorm(vec): It evaluates the norm (4>\(f>) of the state associated with the 
object vec 

[ Input: Matrix, Output: number or expression, Calls: Quantavo, LinearAlgebra] 

StateNormalize(vec/matcol/mat): This procedure returns a normalized 
state such that (ip\ip) = 1 for vec and Tr(p)=l for matcol or mat. If the state was 
not normalized to begin with it will print: "the state is not normalized" and then 
return the normalized state. 

[ Input: Matrix, Output: Matrix and printed answer, Calls: Quantavo, LinearAlgebra, 
IsNormalized, Little Trace] 

StatePartialTranspose(vec/matcol/mat,s): This procedure works only 
for bipartite states (two modes only). If your state is not bipartite use Traceout to 
trace out the other modes. The parameter 's' specifies if the partial transposition must 
be made with respect to the first or to the second mode. The output is the partially 
transposed matcol or mat. 

[ Input: (Matrix, 1 or 2), Output: Matrix, Calls: Quantavo, LinearAlgebra] 

T 

TensorProduct(vec/matcol, List, vec/matcol, List): This will 

return the tensor product between modes [Lis£a] of state A and modes [Lists] of 
state B. It therefore implements the operation, (p, list A, a, listB) — > pu st A ®o"iists- 
Note that the dimension of each state and its corresponding list must be equal. So for 
example we can declare two states, 

Sq:=SqueezedVac(2, 3, £); 

Fock2:=Matrix([[7, [1,1]], [gamma, [2,2]]]); 



7 [1,1] " 
7 [2,2] _ 



1 [0,0] 

e [i.i] 

. e [2,2] . 

and after using 
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TensorProduct (Sq, [1,4] ,Fock2, [2,3] ) ; 



we obtained the tensored vec object, 



7 


[0,1,1,0] - 


7 


[0,2,2,0] 


£7 


[1,1,1,1] 


£7 


[1,2,2,1] 


C 2 7 


[2,1,1,2] 


. C 2 7 


[2,2,2,2] . 



[ Input: Matrix, List, Matrix, List, Output: Matrix, Calls: Quantavo, LinearAlgebra] 

TensorVac(vec/matcol/mat, m): This procedure tensors 'm' vacuum modes 
with an existing state. The new modes are in a product state with the original state. 
If it is a vec object it will therefore do the transformation \ip) — > \tp) ® |0}® m and 
p — > p® (|0){0|)® m for mat and mated. The procedure will also transform the 
global variable K — » K + m. 

[ Input: (Matrix, whole number), Output: Matrix, Calls: Quantavo, LinearAlgebra] 

Traceout(matcol,i): This procedure takes the partial trace of a sparse den- 
sity matrix with respect to mode i. Its input is an object of type matcol 
[ Input: (Matrix, whole number), Output: Matrix, Calls: Quantavo, LinearAlgebra] 

Trim(Vector/vec/matcol): This procedure eliminates all the non zero en- 
tries of a vector of dimension d K and converts it into an object of type vec. It also 
deletes all non-zero entries in objects of type vec and matcol 
[ Input: Vector/Matrix, Output: Matrix, Calls: Quantavo, LinearAlgebra] 

u 

UnitaryEvolution(Unitary Matrix, vec/matcol): Whether we 

have just built a unitary matrix with BuildUnitary or we have an arbitrary unitary 
matrix, we can use this procedure as follows: 



V:= UnitaryEvolution(U, vec/matcol): 



where U is the unitary matrix of dimension K x K and our state is described by 
a vec or a matcol object. This will effectively implement the mode transformation: 

a — Ua 

and return the vec or matcol after the transformation. (Note that for matcol it also 
does b' = U^b ) 

[ Input: (Matrix, Matrix), Output: Matrix, Calls: Quantavo, LinearAlgebra] 

V 

Vac(Nr. of modes): This procedure returns a vacuum vec with the num- 
ber of specified modes. For example Vac(3) will return [ 1 [0,0,0] ] 
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[ Input: positive natural number, Output: Matrix, Calls: Quantavo, LinearAlgebra] 

Vec2ma~t(vec): It converts a vec object into a mat object. The procedures is 
independent of d and K. 

[ Input: 2 column Matrix, Output: Matrix, Calls: Quantavo, LinearAlgebra] 
vec2ma~tC0l(vec) : It converts a vec object into a matcol object effectively doing 

y~]atn\n) — > ^ anafh\n){m\ 

n fifh 

Note that no normalization is implemented. The procedure is independent of d and 
K. 

[ Input: 2 column Matrix, Output: 3 Column Matrix, Calls: Quantavo, LinearAlgebra] 

vec2poly(vec) : It converts a vec object into a poly object 

[ Input: 2 column Matrix, Output: Polynomial, Calls: LinearAlgebra, Quantavo] 

VectorModes(i): this procedure takes as input an integer between 1 and d K 
and outputs a list with the equivalent number of photons in each mode. This way, 
VectorModes(l) = [0,0,0] if it is a 3 mode state and VectorModes(2) = [l,0,0] if K = 3 
and d — 2 for example. (Note that d and K need to be defined) 
[ Input: integer, Output: List, Calls: Quantavo] 

Vec~torRow(Indi, d) : d (maximum number of photons) and K (number of modes) 
need to be specified. This is the inverse procedure of VectorModes. When a list with 
the number of photons in each mode (for instance Indi = [0,1,0]) and the value of d 
are given, it outputs the row number that corresponds to it in a vector type object). 
This way, for d = 2, K = 3 VectorRow([0,l,0],2)=3. 
[ Input: (List,d), Output: whole number, Calls: Quantavo] 

Quantavo also uses the local procedures: 

multiplymatcol, multiplymatcolvec, Tribullesmatcol, Tribullesvec, VectorModes, Vec- 
torRow, indexvec, modesvec, vecBS, myvecBS, matcolBS, mymatcolBS, Projectvecvec, 
Projectmatcol, Dvec, Dmat, Dmatcol, barra,histo; 

B Copyright and Disclaimer 

Copyright (c) 2008 Alvaro Feito Boirac. 

This is the Module QUANTAVO, a toolbox for Quantum Optics calculations that 
can be used in Maple™ (Waterloo Maple Inc.) 

It is released under the GNU General Public License v3 which can be obtained at 
http://www.gnu.org/licenses/gpl.html . Please acknowledge its use if used to es- 
tablish results for a published work. If you make any improvements or find any bugs 
the author will be thankful if you can let him know. 

DISCLAIMER: This Module is distributed in the hope that it will be useful, 
but WITHOUT ANY WARRANTY, without even the implied warranty of MER- 
CHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 

The author can be contacted at abl805@imperial.ac.uk 
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